 function [ysim,stmat]=simgensys(shocks,G1,C,impact,zmat,Y0)
%function [ysim,stmat]=simgensys(shocks,G1,C,impact,zmat,Y0);
% Simulate the path of the variables given the output of gensys.m  
% Inputs 
% ------
% shocks    [NX x T ] vector of shoclk
% G1        Output from gensys 
% C         Constant, Output from gensys
% impact    Output from gensys
% zmat      (optional)  Matrix multiplies the state vector resulting in a matrix of 
%           observables 
%           If not defined, set to the identity matrix 
% Y0        (optional) first state 
%           If not defined, set to zero
%
% Alejandro Justiniano 5/27/03 
% Modified 5/20/04 Added the constant
% Oct 2 2007  Matches what would obtain from KFILTER if fed the 
% simulated shocks here 
% ------------------------------------------------------
if ~isequal(size(C)',[size(G1,1);1] ) ; 
    error('Constant AR matrix dimension mistmatch') 
end 
if nargin < 5; 
    zmat=eye(size(G1,1));
    disp('Simulate all variables in the system') ; 
end
T=size(shocks,2); 
ysim=zeros(size(zmat,1) ,T ); 
stmat=zeros(size(G1,1),T); 
if nargin < 6; 
    stmat(:,1)=zeros(size(G1,1),1)+impact*shocks(:,1); 
else 
    stmat(:,1)=Y0;
end
ysim(:,1)=zmat*(stmat(:,1)+C); 
%disp(['Begin simulation for ',num2str(T),' periods '] ); 
for ii=2:T ;     
    stmat(:,ii)= G1*stmat(:,ii-1)+impact*shocks(:,ii);
    ysim(:,ii)=zmat*(stmat(:,ii)+C);  
end 
